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THE AXIALLY SYMMETRIC POTENTIAL FLOW ABOUT 
ELONGATED. BODIES OF REVOLUTION 


by 


L. Landweber 


ABSTRACT 


An iteration formula for Fredholm integral equations of the first kind is ap- 
plied in two new methods for obtaining the steady, irrotational, axisymmetric flow of 
an inviscid, incompressible fluid about a body of revolution. In the first method a 
continuous, axial distribution of doublets is sought as a solution of an integral equa- 
tion of the first kind. A method of determining the end points and the initial trends 
of the distribution, and a first approximation to a solution of the integral equation are 
given. This approximation is then used to obtain a sequence of successive approxima- 
tions whose successive differences furnish a geometric measure of the accuracy of an 
approximation. When a doublet distribution has been assumed, the velocity and pres- 
sure can be computed by means of formulas which are also given. 

In the second method the velocity is given directly as the solution of an inte- 
gral equation of the first kind. Here also a first approximation is derived and applied 
to obtain a sequence of successive approximations. In contrast with the first method, 
which, in general, can give only an approximate solution, the integral equation of the 
second method has an exact solution. 

Both methods are illustrated in detail by an example. The results are com- 


pared with those obtained by other well-known methods. 


INTRODUCTION 


HISTORY > 


The determination of the flow about elongated bodies of revolution 
is of great practical and theoretical importance in aero- and hydrodynamics. 
Such knowledge is required in connection with bodies such as airships, tor- 
pedoes, projectiles, airplane fuselages, pitot tubes, etc. Since it is well 
known that for a streamlined body, moving in the direction of the axis of sym- 
metry, the actual flow is very closely approximated by the potential (inviscid) 
flow about the body,? numerous attempts have been made to find a convenient 
theoretical method for obtaining numerical solutions of the potential flow 
problem. 

At first the problem was attacked by indirect means. In 1871 
Rankine* showed how one could obtain families of bodies of revolution of known 
potential flow, generated by placing several point sources and sinks of vari- 
ous strengths on the axis. This method was extended and used by D.W. Taylor® 
in 1894 and by G. Fuhrmann* in 1911. The latter also constructed models of 
the computed forms and showed that the measured distributions of the pressures 
over them agreed very well with the computed values except for a small region 
at the downstream ends. More recently, in 1944, the Rankine method was em- 
ployed by Munzer and Reichardt’ to obtain bodies with flat pressure distribu- 
tion curves, and a further refinement of the technique was published by 
Riegels and Brandt.” Most recently the indirect method has been employed to 
obtain bodies generated by axisymmetric source-sink distributions on circum- 
ferences, rings, disks,and cylinders. This development, which enabled bodies 
with much blunter noses to be generated, was initiated by Weinstein | in 1948 
and continued by van Tuyl > and by Sadowsky and Sternberg” in 1950. 

A method of solving the direct problem, i.e., to determine the flow 
over a given body of revolution, appears to have been first published by 
von Kdrmdn?® in 1927. von Karman reduced the problem to that of solving a 
Fredholm integral equation of the first kind for the axial source-sink distri- 
bution which would generate the given body, and solved the integral equation 
approximately by replacing it with a set of simultaneous linear equations. 
Although this method has limited accuracy and becomes very laborious when, for 
greater refinement, a large number of linear equations are employed, neverthe- 
less it is the best known and most frequently used of the direct methods. A 

‘modification of the von Kdrmdn method was published by Wijngaarden?? in 1948. 


*References are listed on page 59. 


An interesting attempt to solve the direct problem was made by 
Weinig?® in 1928. He also formulated the problem in terms of an integral equa- 
tion for an axial doublet distribution which would generate the given body, 
and employed an iteration formula to obtain successive approximations. Since 
the successive approximations diverged, the recommended procedure was to extra- 
polate one step backwards to obtain a solution. 

In 1935 an entirely different approach, in which a solution for the 
velocity potential was assumed in the form of an infinite linear sum of orthog- 
onal functions, was made by Kaplan?* and independently by Smith.+* The coeffi- 
cients of this series are given as the solution of a set of linear equations, 
infinite in number. In practice a finite number of these equations is solved 
for a finite number of coefficients, and Kaplan has shown that the approximate 
solution thus obtained is that due to an axial source-sink distribution which 
is also determined. A simplification of Kaplan's method by means of addition- 
al approximations was proposed by Young and Owen?® in 1943. 

It appears to be generally agreed, by those who have tried them, 
that the aforementioned methods are both laborious and approximate. Thus, ac- 
cording to Young and Owen:+° 

"In every case, however, the methods proposed are laborious 

to apply, and the labour and heaviness of the computations 

increase rapidly with the rigour and accuracy of the proc- 

ess. Inevitably, a compromise is necessary between the 

accuracy aimed at and the difficulties of computation. All 

the methods reduce, ultimately, to finding in one way or 

another the equivalent sink-source distribution, and it is 

this part of the process which in general involves the 

heaviest computing." 

Furthermore, a fundamental objection is that only a special class of bodies of 
revolution can be represented by a distribution of sources and Sinks on the 
axis of symmetry. According to von Karman;1?° 

"This (representability by an axial source-sink distribu- 

tion) is possible only in the exceptional case when the 

analytical continuation of the potential function, free 

from singularities in the space outside the body, can be 

extended to the axis of symmetry without encountering 

Singular spots." 

The dissatisfaction with these methods is reflected by the continuing attempts 
to devise other procedures. 

A new method published by Kaplan*® in 1943 is free of the assumption 
of axial singularities and appears to be exact in the sense that the solution 
ean be made as accurate as desired, but the labor required for the same ac- 
curacy appears to be much greater than by other methods. The application of 
the method requires that first the conformal transformation which transforms 


the given meridian profile into a circle be determined. The velocity poten- 
tial is then expressed as an infinite series whose terms are universal func- 
tions involving the coefficients of the conformal transformation. Kaplan?® 
has derived only the first three of these universal functions. 

Cummins of the David Taylor Model Basin is developing a method based 
on a distribution of sources and sinks on the surface of the given body. This 
method is also exact, but the labor involved in its application has not yet 
been evaluated. 

Another exact method, based on a distribution of vorticity over the 
surface of the body, is being developed by Dr. Vandry of the Admiralty Re- 
search Laboratory, Teddington, England. The methods of both Cummins and 
Vandry lead to Fredholm integral equations of the second kind, which can be 
solved by iteration. 

The present writer has developed two new methods, an approximate one 
in which an axial doublet distribution is assumed, and an exact one based on 
a general application of Green's theorem of potential theory. Both methods 
lead to Fredholm integral equations of the first kind for which a solution by 
iteration has been discussed by the author.?’ Indeed, the consideration of 
this iteration formula was initiated in an attempt to find more satisfactory 
solutions of the integral equations of von Karman?® and Weinig.22 These new 
methods will be presented, and, by application to a particular body, compared 
with other methods from the point of view of accuracy and convenience of 
application. 


FORMULATION OF THE PROBLEM 


We will consider the steady, irrotational, axially symmetric flow 
of a perfect incompressible fluid about a body of revolution. Take the x-axis 
as the axis of symmetry and let x, y be the coordinates in a meridian plane. 
Denote the equation of the body profile by 


3p = 38((33) [1] 


Since the flow is irrotational there exists a velocity potential ¢ 
which, for axisymmetric flows, depends only on the cylindrical coordinates 
x, y and satisfies Laplace's equation in cylindrical coordinates 


dx (Yee) + ay (98g) 0 (21 


Also, since the flow is axisymmetric, there exists a Stokes stream function 
y(x, y) which is related to the velocity potential by the equations 


f) 3 f) 3 
aes ay Se oe gee [3] 


It is seen that Equation [2] may be interpreted as the necessary and suffi- 
cient condition insuring the existence of the function w. As is well known, 
w is constant along a streamline and, considering the surface of revolution 
generated by rotation of a streamline about the axis of symmetry, 2my may be 
considered as the flux bounded by this surface. On the surface of the given 
body and along the axis of symmetry outside the body we have w= 0. wy satis- 
fies the equation 


which is obtained by eliminating ¢ between Equations [3]. 

The velocity will be taken as the negative gradient of the velocity 
potential. Let u, v be the velocity components in the x, y directions, Then 
by [3], we have 


Ep Op el walked. 
OOK == [5] 
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For a uniform flow of velocity U parallel to the x-axis we have 
@= -Ux, p= -5 Uy® [7] 


The boundary condition for the body to be a stream surface may be 
written in various ways. If the body is stationary the boundary condition is 


OD (2, VY =O [8a ] 
or, equivalently, 


(38), - 0 [8b] 


where the derivative in [8b] is evaluated on the surface of the body in the 
direction of the outward normal to the body. If the body is moving with veloc- 
ity V parallel to the x-axis the boundary condition becomes 


(eel = -V cos 8 [9] 


where B is the angle between the outward normal to the body and the x-axis. 
It is desired to obtain a solution of [2] or [4] which satisfies the 
boundary conditions [7] at infinity and [8] or [9] on the body. 


METHOD OF AXIAL DISTRIBUTIONS 


SOURCES AND SINKS 


The potential and stream functions for a point source of strength M 


situated on the x-axis at x = t are 
M NM x-t 
0 y= M(-1 + 25) [10] 


where 


Tey at) aaet aya (a) 


If the sources are distributed piecewise-continuously along the x- 
axis between the points a and b (see Figure 1) with a strength w(x) per unit 
length, the potential and stream functions are 


Figure 1 - The Meridian Plane 


As is well known, Rankine bodies are obtained by superposition of 
these flows with a uniform stream so as to obtain a dividing streamline begin- 
ning at a stagnation point. Without loss of generality we may suppose this 
uniform stream to be of unit magnitude. This dividing streamline is the pro- 
file of the Rankine body for which, by [7], the stream function is 


abe pier ay) ees 
b= -ty + fate) (-1 + =*) at [14] 
The boundary condition, Equation [8a], then gives as the implicit equation for 
the body 


fae) (21 +%)at=gy? [15] 


where now y* = f(x) and r® = (x-t)® + f(x). In order to obtain a closed body 
the total strength of sources and sinks must be zero, i.e., 


[med se 0 


a 


In this case [15] becomes 


[wey 28 at = 3 [15a] 

In general [15a] cannot be solved explicitly for f(x) when u(t) is 
given. A practical procedure for obtaining f(x) for a given x is to evaluate 
the integral numerically for various assumed values of f(x) and to determine 
the value which satisfies [15a] by graphical means. 

When f(x) is prescribed [15a] may be considered as a Fredholm inte- 
gral equation of the first kind for determining the unknown function u(t). 
This equation will not be treated. Indeed it will be shown that, when con- 
tinuous distributions are considered, it is a special case of the more general 
equation for doublet distributions which will now be derived. 


DOUBLET DISTRIBUTIONS 


Let m(x) be the strength per unit length of a continuous distribu- 
tion of doublets along the x-axis between the points a and b (see Figure 1). 
The potential and stream functions may be taken as 


ae [ m(t) = at [16] 


a TP 


and 


is Al, mE) at 017] 


The stream function for a Rankine flow now becomes 


b 
y= bye ey? [| MES ot [18] 


r 


Hence the boundary condition, Equation [8a], gives 


PInGb ames suc 
fae at -5 [19] 


Here again Equation [19] may be considered as an implicit equation for the 
Rankine body when m(t) is given, or as a Fredholm integral equation of the 
first kind when the body profile y* = f(x) is prescribed. 

In order to show the relation between the source and doublet distri- 
butions in Equations [15a] and [19], integrate by parts in [19]. We have 


b 2 6 6 
SE ae) o-x dm x-t 
[mer 3G at = mt) © +48 = at 
Hence [19] may be written as 
m(t ) fe [ ° dm x-t at oils 2 [20] 
Teas te 5 Gig 2 


The interpretation of Equation [20] is that a doublet distribution of strength 
m is equivalent to a source-sink distribution of strength dm/dt together with 
point sources of strength m(a) and -m(b) at the end points. Hence source-sink 
distributions are completely equivalent only tc those doublet distributions 
which vanish at the end points. This justifies the remark in the previous 
section that the integral equation for the doublet distributions is more gen- 
eral than that for the source-sink distributions. 


MUNK'S APPROXIMATE DISTRIBUTION 


Munk?® has given an approximate solution of Equation [19] for elon- 
gated bodies. His formula may be derived as follows. At a great distance 
from the ends of a very elongated body, the integrand of [19], m(t)/r°, will 


peak sharply in the neighborhood of t = x. In the range of the peak, in which 
the value of the integral is principally determined, m(t) will vary little 
from m(x). Also, only a small error will be introduced by replacing the lim- 
its of integration by -~ and+co. Hence, as a first approximation to a solu- 
tion of [19], try 


mix)f = 3 [21] 


We obtain 


n, (x) = 7 ¥ [22] 


a distribution proportional to the section-area curve of the body. This ap- 
“proximation was independently derived by Weinig?? who employed it as the first 
step in a divergent iteration procedure. It has also been rediscovered by 
Young and Owen?® and Laitone’® who have shown the accuracy of the approxima- 
tion for elongated bodies by several examples. 

It is apparent from its derivation that [22] also gives the asymptot- 
ic radius of the half-body generated by a constant axial dipole distribution 
extending from a point on the axis to infinity. It is readily seen that this 
distribution is equivalent to a point source at the initial point. 

As a refinement to Munk's formula, Weinblum®® has used the approxi- 
mation 


m, (x) =, Cy* [23 ] 
where C is a factor obtained by comparison of the distributions and section- 
area curves of several bodies. Weinblun's factor bears an interesting rela- 

A 
tion to the virtual mass of the body. This is seen by considering the expres- 


sion for the virtual mass k 4 in terms of the mass of the displaced fluid A 
and the totality of the AoUBIStE, fc nebe = PBS DES) 


b 
k,A = 4mp mdx - A [24 ] 


where k, is designated the longitudinal virtual mass coefficient, and p is the 
density of the fluid. But, from [23], 


b b 
4p | m, ax = po | my2dx = 4CA 
a a 
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since, for elongated bodies, a and b very nearly coincide with the body ends. 


Hence 


c = +k) (25 ] 


In practice an approximate value of k, may be taken as that of the 
prolate spheroid having the same length-diameter ratio as the given body. The 
values of k, for a prolate spheroid may be computed from the formula®* 


4 xin (a + VaF-1) - VAF-1 [26] 
2 VP =T = aln (A+ VA? =1 


where A is the length-diameter ratio. Hence 


1 3/2 
oie a [27] 
2 Vr2-1 - Aln (a + Vie 


The values of k, versus d have also been tabulated by Lamb” and graphed by 


Munk .° 


END POINTS OF A DISTRIBUTION 


A difficulty in determining the doublet distribution from Equation 
[19] is that the limits of integration, a and b, are also unknown. In the 
method of von Karman?® the end points are arbitrarily chosen; Kaplan?® takes 
the end point of the distribution midway between the end of the body and the 
center of curvature at that end. 

Kaplan based his choice on a consideration of the prolate spheroid. 
Thus the equation of the spheroid of unit length and length-diameter ratio i, 
extending strom: xs == OM GOlexd illus 


1 
y° = (x - se) [28] 
x 
The radius of curvature at x = 0 is then +5. The exact doublet distribution, 
however, extends between the foci of the spheroid which are situated at dis- 


tances 
NING ee) 
OX 


from the end points. Hence the error in Kaplan's assumption, 


NN ge Silage el ( EINE 
2a Xe alone 2r? 


diminishes rapidly with increasing 2. 

For the half-body generated by a constant doublet distribution (a 
point source), Kaplan's assumption gives a poor approximation. Let a® be the 
strength of the distribution. Then it can easily be shown from [19] that the 
source is at a distance a from the end of the body (stagnation point), and 
that, if the origin is chosen at the latter point, the equation of the half- 
body is 

2 2 3 
(Sar a) tae a) (291 
Hence the radius of curvature at the end 1s 52, so that Kaplan's Seas ee 
for the start of the distribution gives za. This is in error by 32. 
An approximate method for determining the end points of a distribu- 


tion and its trends at the ends is given in Appendix 1. The given profile is 
assumed to extend from x = 0 to x = 1 and to have the equation 


= 2 3 
Yr 2 Ge BLS cb AES A oon [30] 


The doublet distribution is assumed to extend from x = a to x = b, So that 
O<a<b<1, with a near O and b near 1, and to have the equation 


m(x) = ¢, + c,x + Ce fas ce [31] 


0 
Only the trends of the distribution near the origin are discussed in Appendix 
1. It is clear, however, that by means of a linear transformation the equa- 
tion of the given profile can be expressed so that the end points of the body 
exchange their roles. Hence the results in Appendix 1 can be applied to either 
end of the body. 

The method of Appendix 1 consists essentially of expanding the inte- 
gral in [19] about the origin and equating powers of t on the two sides of the 
equation to obtain a series of equations in the unknowns a, Cy» Cis Cc 


2 


2 
By applying the first four of these equations an approximate solution is ob- 


tained in the form 


ied} 


a=— [32] 


where 


and @ 


where 


12 


“Ha? [309 - 37a° + 1200 - 96 + Zila, + 2haGar - 15a + 16 - Ha)| al 
a [150° - 1680° + 5120 - 384 + 96a, + 48a(50° - 24a t 24 - 6a) [34] 


2 


o,D = -4[(a - 4)2(a - 1) + 4a,] [35] 


D = 2(9a° - 94a? + 2720 - 192) + 8[(a sh2ia 2 a) 4 4a] in a 
r See > Galen So Aehi’= «> QIMGH > Tcts))) > See, [36] 
- 96a? (5a? - 24a + 24) + 516a°a, 

is a root of the seventh-degree polynomial 


2 2 2 
A + a,B + a0 + aia,D + a,k + a,a,F + aa,G + a,aiH = (0) [37] 


A(a) = ala - 4)? (5a* - 830° + 2880 - 368a + 128) 
Bla) = 72(a - 4)? (50° - 250° + 40a - 16) 

Cla) = Yala - 4) (530° - 148a + 128) 

D(a) = -288(a@ - 4)(5a* - 160 + 16) 

E(a) = -96a(3e - 4) 

F(a) = 1152(2e - 3) 


48a(3a - 8) 


Q 
R 
i 


H(a) = -1152(@ - 3) 


The solution gives, for the initial doublet strength at x = a, 


m(a) 


[la - 4) (a? - 120 + 16) + 4Bala - ¥)(a - 2) + 16a, - 96aa,)] [39] 


When ay, aj, 4, .-- are all small in comparison with unity, an ap- 


proximate solution for a@ is 


W 


/ 

b aut S 

Qs Web a =c ey le a 2a [40 ] 
a-4+a, ifa <0 [41] 


and, to the same order of approximation, 


a a a az 
ax) =O +5 +52 ne \(-e ty) [42] 
and 
a a a 
mal) -3(1 + = i = in a VECEo | ue ey, 20 [43 ] 
(ay sO, i a, < © [44 ] 


It is seen that Kaplan's assumption that a = 4 gives the principal 
term of the solution in [40] or [41]. The form [42] immediately suggests a 
modification and refinement of the Munk-Weinblum approximation, Equation [23], 
which will be considered in the next section. 

A graphical procedure for finding the roots a@ of Equation [35] is 
also given in the Appendix. For this purpose the functions A(a), Bla),... Hla) 
are tabulated in Table 10. 


AN IMPROVED FIRST APPROXIMATION 


According to its derivation the Munk approximation could be expected 
to be useful only at a distance from the end points of a distribution. It was 
seen, however, Equation [42], that under certain circumstances a distribution 
which was a suitable approximation for the nose and tail of a body also ap- 
peared as a generalization of the Munk-Weinblum approximation, [23]. This 
suggests a procedure for obtaining an improved approximate distribution. 

It is desired to obtain a distribution m(x) which satisfies the fol- 
lowing conditions: 


(a) m(x) assumes known values m, and m, at the distribution limits a and 
be pent; 


m(a) = m m(b) = M,, [45 ] 


(b) m(x) is nearly equivalent to the Munk-Weinblum approximation [23] at 
a distance from the distribution limits, i.e., 


may 1Cvs for a <x <<)b 
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(c) m(x) satisfies the virtual mass relation [24] which may be written 
in the convenient form 


[ m(xax = 7 + i) { yPax (46 ] 
a 0 


It is readily verified that Condition (a) is satisfied by the 
distribution 


m(x) = Cy° +e + ex [47] 
where 
©) = pea (DM, - om, + Caf, - of,)| [48 ] 
and 
e = pal" -m, + C(f, - f,)] [49] 


If the linear term e, + ex in [47] is small in comparison with m(x) at a dis- 
tance from the ends, then Condition (b) is also satisfied. Finally, Condition 
(c) can be satisfied by a proper choice of C in [47]. This is accomplished by 
writing m(x) in the form 


m(x) = Cy? - px es = f,) py ee Be 
substituting it into Equation [46], and solving for C. We obtain 


aL aie i 2 dx - 1(b-a)(m +m, ) 
mM dhs 2 ab 
b 2 ae : 
IL y>dx - o(b-a)(f,+f,) 


[50] 


SOLUTION OF INTEGRAL EQUATION BY ITERATION 


Now that we have derived a good first approximation to the doublet 
distribution function in the integral equation [19], it would be very desir- 
able to apply it to obtain a second, closer approximation. This can be accom- 
plished by means of the iteration formula which we will now derive. 

Let m,(x) be a known first approximation and yp, (x) the corresponding 
values of the stream function won the given profile y> = f(x). Then, from 
Equation [18], 


1 


(t) 
Y,(x) = - f(x) + r(x) [ - dt 


[51] 


Thus wy, (x) is a measure of the error when m, (t) is tried as a solution of the 
integral equation [19]. If m(t) is a solution of [19], Equation [51] may be 
written in the form 


b m, (t)-m(t) 


3 dt [52] 


p(x) = (x) | 


a 10 


But, on the same assumptions as were used to derive Munk's approximate distri- 
bution, Equation [22], we obtain as an approximate solution of the integral 
equation [52] 


m, (x) - m(x) = 3 (x) [53] 


or, denoting the new approximation to m(x) by m, (x), 


m,(x) = m,(x) - 3p, (x) [54] 
Hence, from [51] 
om, (t) 
m(x) = m, (x) deez =|, ae | [55] 
ces 


Since the foregoing procedure can be repeated successively, we obtain the iter- 


ation formula 


1 : om, (t) 
ms 44 (x) = m, (x) + + 1003 -| 2 dt [56] 
and 
m,,4(*) - m,(x) = -5 py, (x) [57] 


It is seen that Vs is the value of the stream function on the given 


profile corresponding to the th approximation m,(x) and hence serves as a 


measure of the error when m,(t) is tried as a Sonesta of the integral equa- 
tion [19]. 

Although successive approximations to m(x) may be computed directly 
from [56], an alternative form, which is both more convenient and more signif- 


icant, will now be derived. From [56] we may write. 
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m, (x) = m,_,(x) + + t(x) Pa | [56a J 
Hence, deducting [56a] from [56] and making use of [57], we get 
B(x) = Wy yx) see at [58] 
Also, from [57] we obtain 


i 
my (x) = m (x) - > Sy, (x) (591 


Thus, in order to obtain ms ,4(x), we first assume an m, (x), then determine 
, (x) from [51]. p(x), p(X), ... can then be successively obtained from 
[58], and finally m, 1 (X) from [59]. 

It has been stated that the magnitude of yp, (x) is a measure of the 
proximity of m, (x). This property of W, (x) can be given a geometrical in- 
terpretation. Corresponding to the distribution m,(x) there is an exact 
stream surface on which the stream function B, (x, a) = Oy leet An, be the 
distance from a point (x, y) on the given body to this exact stream surface, 
measured along the normal to the given body, positive outwards. Let uy be 
the tangential component of the flow along the body. Then we have 


1 Ap, (x,y) 
Poewes a) 8! Onan ati 


But Ay = -y, (x), since W(x, y) = 0 on the exact stream surface. Hence 


x 

ies 2 : [60] 
Ss 

Since, for an elongated body, UES 1, except in the neighborhood of the stag- 
nation points, it is seen that , (x) enables a rapid estimate to be made of 
the variation from the desired profile of the exact stream surface correspond- 
ing to m, (x). This is an important property because it can be used to monitor 
the successive approximations. Thus, the sequence Wp; (x) can be terminated 
when An, becomes uniformly less than some specified tolerance; or, since there 
is no assurance that the infinite sequence p, (x) converges, the sequence can 
conceivably give useful results even without convergence if it is continued as 


7 


long as An, decreases on the average, and is terminated when the error begins 
to increase and grows to an unacceptable magnitude at some point along the 
body. The strong similarity between these remarks and the discussion follow- 
ing Theorem 2 of Reference 1/ should be noted. 

There is also a strong similarity between the iteration formula of 
Reference 17 whose convergence was thoroughly discussed, and the present equa- 
tion [56]. An essential difference between the iteration formulas is that the 
former employs the iterated kernel of the integral equation, the latter does 
not, so that the convergence theorems of Reference 1/ are not applicable. Nev- 
ertheless, it is proposed to use the form in [56] (or the equivalent iteration 
formula [58]), for the following reasons: 


a. The labor of numerical calculations would be greatly increased by 
iterating the kernel, and even then only convergence in the mean would be 
guaranteed (Theorem 4 of Reference 17). 


b. The physical derivation of Equation [56] indicates that at least the 
first few approximations should be successively improving. 


ce. The successive approximations are monitored so that the sequence can 
be stopped when the error is as small as desired or, in the case of initial 
convergence and then divergence, when the errors begin to grow. 


VELOCITY AND PRESSURE DISTRIBUTION ON THE SURFACE 


When an approximate doublet distribution m, (x) has been obtained, 
the velocity components u, v can be computed from the corresponding stream 
function [18] 


bm, (t) 
y, (x,y) = y* f as dt - 1] [61] 


from which, in accordance with Equations [5] and [6], 


and 
-X 


t 
v = 3y ee m, (t)dt [63 ] 
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On the given surface we have, from [61], 


(t (x) 
[A een ee [64] 
Dig y* (x) 
where now 
> ce (ee) + ses) [65 ] 
Differentiating [64] with respect to x gives 
Wi(x) 2, (x)y"(x) 
jf Se nei oS [66] 
a or y> (x) y? (x) 
Hence, from [62] and [64] we obtain 
> m, (t) 2p, (x) 
= dt - [67] 
i I r> f(x) 
and, from [63], [66], and [67], 
see) 
v = uy!(x) + oes [68] 
y(x) 


where the primes denote differentiation with respect to x. Equations [67 ] 
and [68] are the desired expressions for u and v. If the approximation m,(t) 
is very good, the contributions of the error function w,(x) should be very 
small. It is interesting to note that the form of Equation [68] shows the 
deviation of the resultant velocity from the tangent to the given body. 

Bernoulli's equation for steady, incompressible, irrotational flow 
with zero pressure at infinity now gives the pressure distribution p, 


= 1 - (u? + v?) [69] 


ald 


where q is the stagnation pressure. 
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NUMERICAL EVALUATION OF INTEGRALS 


In order to perform the iterations in Equations [56] and [58] and to 
compute the velocity distribution it will frequently be necessary to evaluate 
integrals of the form 


3 


ti) b 
mt) gt ana { Se at 
a op B ip> 


where 


Po (c=) ieee (x) 


Because in this form these integrals peak sharply in the neighborhood of t = x, 
especially when the body is elongated, they are consequently unsuited for nu- 
merical evaluation. 

A more suitable form can be obtained by means of the following trans- 
formation. Let (x, y) be the coordinates of a point on the body, t the ab- 
scissa of a point on the axis, 6 the angle between a line joining these two 
points and the x-axis; see Figure 1. Then 


x - t = y(x) cot 6 [70] 


We may now transform the integrals so that 6 becomes the variable of integra- 


tion. Then 


b y? B 
{z nteere = | ae) ened [71] 
an a 
and 
b y4 B 
== m(t) = i m(t) sin® @ d@ [72] 
G12 a 
where 
= Ni 2 ae 
a = arctan >=, = arctan >> [73] 


An alternate procedure, which eliminates the peak without a trans- 
formation of variables, is the following. We have 


(Ss m(t)at = [ S[pt)-a00] aes) (3 a 
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and 


ie m(t)dt = [ %[n(er-moo] et + m(x) "Se at 
a r> la rs L r> 


Hence 


ls m(t)dt = [5 [mer -noo) at + m(x)(cos@ - cos) [71a] 


if SE ae ee = [ E [mct)-moxr]at 


r> rp? 
+ m(x) [eos ao -cosB = (costa - cos?) | [72a ] 


Gauss! quadrature formula is a convenient and accurate method of 
evaluating these integrals. The formula may be expressed in the form 


1 
[Pleas - & RaF(Ena) [74] 


where the gy are the zeros of Legendre!'s polynomial of degree n and the Rai 
are weighting factors. These have been tabulated?® for values of n from 1 to 
16. These numbers have the properties 


R eral Ga OG aa sing [75] 


ni > ®y nit 
The value of the integral given by Formula [74] is the same as could be ob- 
tained by fitting a polynomial of degree 2n-1 to F(x). The values of Re and 
are tabulated in Table 1 for n = 7, 11, and 16. 

When the limits of integration are w and #, as in Equations [71] 
and [72], Gauss' formula becomes 


Sai 


B B-a ¥ 
[ Fle)ae rene, RagF(8,) [76] 
where 
B- +B 
Oe Ga oe [77] 
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TABLE 1 


ABSCISSAE AND WEIGHTING FACTORS FOR GAUSS! QUADRATURE FORMULA 


-0.949108 0.129485 -0.978229 0.055669 -0.989401 0.027152 
741531 279705 887063 .125580 -QUU575 .062254 
-0.405845 .381830 130152 . 186290 .865631 .095159 
0.417959 .519096 .233194 -755404 124629 

-0.269543 .262805 .617876 .149596 

0 0.272925 -458017 .169157 

R, =R .281604 . 182603 
-0.095013 0.189451 


Gry SSG aeaihy|| May = Aetaey 


$1 > Sn -i41 


al n-it+] 


ILLUSTRATIVE EXAMPLE 


The foregoing considerations will now be applied to a body of rev- 
olution whose meridian profile is given, for -1< x S11, by, 


y? = f(x) = 0.04(1 - x*) [78] 


The body is symmetric fore and aft, has a length-diameter ratio A = 5, and a 
prismatic coefficient 


1 
ae | (1 - x*)ax = 0.80 [79] 
0 
By applying to [78] the transformation 
R= 26S. yoaey [80 ] 


We obtain the equation for the geometrically similar body of unit length, for 
0oS¢@ S14, 


n? = 0.08(€ - 362 + 4é® - 26%) = 0.086(1 - €)(26? - 2€ + 1) [81 ] 


Ze 


We will also need the slope of the profile which, from [78], is 


fell (5s) ees Oe [82 ] 


Die 
J 2y (1-x4)272 


The profile and f(x) are graphed in Figure 2. 


Figure 2 - Graphs of y(x) and y*(x) for y?(x) = 0.04(1 - x‘) 


First let us find the end points of the distribution. We have, from 
isl, 2, = Oo, A, 2S Ooty A, = 0.32. The approximate formula [40] then 
gives a = 3.68 or 3.84, whence a = > = 0.0217 or 0.0208. An examination of 
the complete polynomial [37] with the aid of Table 10 shows that its zeros oc- 
Ge Be OS 35655, Boel), Wale iba Wne application of Table 10 to determine these 
roots the regions of possible zeros should be determined by inspection, the 
values of the polynomial in these regions calculated from Equation [37] and 
Table 10, and then graphed to obtain the zeros. It is seen that in the pres- 
ent case the approximate formula [40] would have been sufficiently accurate 
for the determination of the roots near a = 4. The solution of the complete 
polynomial equation will always yield an additional large root, corresponding 
to the large root of Equation [131] of the Appendix; in general, however, this 
root should be rejected since as will be shown, the initial doublet distribu- 
tion corresponding to it is less simple than for the roots neare = 4. 


(45) 


The initial behavior of the distributions corresponding to each of 
the three roots, as determined from Equations [33] through [36], and [39], is 
shown in Table 2. It is seen from the table that the distribution for @ = 12.1 
begins with practically a zero value for m(a), with a small negative slope and 
with upward curvature. Since the distribution curve cannot continue very far 
with upward curvature, there must be an inflection point nearby. In contrast, 
the distribution corresponding to the other two roots begin with positive 
slopes and downward curvatures and hence must be considered simpler. Further- 
more, the distribution for the first root is considered simpler than for the 
second since the distribution curves are practically identical except that, 
for the second root, the curve is extended a distance Aa = 0.0011, in the 
course of which m(a) changes from a positive to almost a numerically equal 
negative value. If we take the point of view that the positive and negative 
values of this extension counterbalance each other, the curve without the 
extension, corresponding to the first root, must be considered the simplest. 


TABLE 2 


Characteristics of Initial Distribution 


Hence, for the purpose of obtaining a first approximation, we will 
assume @ = 3.65 and, correspondingly, a = 0.022, m(a) = 0.000022. Often, as 
in this case, the labor of obtaining a and m(a) can be considerably reduced 
by using the less exact equations [40] through [44] instead of [37] through 
[39]. Since, as will be seen, the iteration formulas rapidly improve upon the 
first approximation, great effort should not be expended to determine an 
initial value for m(a). 

The values a = 0.022 and m(a) = 0.000022 have been derived for the 
profile in the é, n-plane. The corresponding values in the x, y-plane are 
a = -0.956 and m, = 0.000088. By symmetry we also have b = -a, m, = m,. 

A first approximation can now be obtained from [47], [48], [49], and 
(50, Stinoa A S 2,9, we have k, = 0.059. Also, from [78]: it = 0.00659, 
[yrax = 0.0640, [ yPax = 0.0637. Hence from [50], C = 0.328. Then, from [48], 
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@ =m, - Cf, = -0.00207; from [49], e, = 0. Finally we obtain from (47 ] 


m,(x) = 0.328y* - 0.00207 [83] 


We can now apply Equation [51] and the iteration formula [58] to ob- 
tain the sequence of functions p, (x). Let us suppose that it is desired to 
obtain a distribution m, (x) whose exact stream surface deviates from the given 
surface by less than one percent of the maximum radius, i.e., An < 0.002. 
Then, by [60], the sequence ; (x) should be continued until Wp, (x) < 0.002 Vf(x) 
for aX x <b, unless the error, as represented by y(x), begins to grow before 
the desired degree of approximation is attained. In the latter case the best 
approximation attainable would fall short of the specified accuracy. 

The integrations in [50] and [51] may be carried out in the form 
[71] in terms of @ defined in [70]. For a fixed (x, y) on the given profile, 
@ and # are first computed from [73]. Then, to apply Gauss! quadrature formu- 
la [76], the interval is subdivided at the points a; given by [77] and the 
integrands evaluated at these points. The corresponding values of t at which 
m, (t) in [51] or p,_(t) in [58] is to be read are, from [70], 


t; = Gn ye COG 9; [70a ] 
Since the values t. and sin @. are used repeatedly in the successive itera- 
tions at a given (x, y), these should be stored in a form convenient for 
application. 

The calculations for obtaining the integration limits @ and @ for 
several values of x are given in Table 3. The values of 6, from [7/7], and the 
corresponding values of R, sin 6. for application of the Gauss 11 ordinate 
formula, and the values of t. from [70a] for each x are entered as the first 
three columns in Tables 5a through 5h, in which are given the calculations for 
(x). 

In order to compute W(x), m,(t) is computed from [83], then 
m,R Sin @ is obtained. These are tabulated in Table 5. Gauss! formula then 
gives fn, singdé. p, (x) is then obtained from [51]; its graph is given in 
Figure 3. It is important to note that m, (t) is obtained by calculation, 
rather than graphically, in this operation. This procedure is recommended 
since it gives greater accuracy in a critical step. In the subsequent opera- 
tions on the y's considerably less percentage accuracy is required, since the 
W's are of the nature of first differences between the m's, so that graphical 
operations are sufficiently accurate. 
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As a check on the accuracy of the integration, wv, (0) was also eval- 
uated by two other means, with the following results: 


From Gauss 7-ordinate formula w,(0) = 0.001258 
From Gauss 1l-ordinate formula vw, (0) = 0.001243 
From exact integration y, (0) = 0.001243 


It is seen that the 7-ordinate formula introduces an error in the fifth deci- 
mal place. 

The first step in the determination of (x) is to read the values 
of w(t) from the graph, Figure 3. wR Sin 6 and Je sin 6d @are then ob- 
tained. w(x) is then given by [58] and graphed in Figure 3. Repeated 
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Figure 3 - Comparison of Error Functions v(x) from 
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Iteration Formula and von Karm4n Method 
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application of this procedure gives w(x) and w(x) which are also graphed in 
Figure 3. The sequence is stopped at p(x) since v, has increased appreciably 
over yp, at x = -0.956. 

Hence, from [59], we have the approximate distribution 


m, (x) = m,(x) - 3 |, (x) + (x) + v, (x)] [84 ] 


to which p(x) is the corresponding error function. The distance An between 
the stream surface for m, (x) and the given profile is seen to be very small; 
the largest error, v, = -0.00007 at x = -0.956, gives a An of about one per- 
cent of the maximum ordinate. A graph of m, (x) is given in Figure 4. For the 
sake of comparison the curves for m, (x) and the original Munk approximation 
rae) are also shown. 


m,(x), M(x), and y7/4 


-10 -09 -08 -07 -06 -05 -04 -03 -0.2 -0.1 Q 


Figure 4 - Comparison of Doublet Distributions m, (x), m (x), 
and Munk!'s Approximation 
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Table 4 shows the calculations for obtaining the velocity components 
u, v from [67] and [68], in which the integrals have been evaluated in terms 
of the polar angle @, according to Equations FAN, (721, eae [Woe kere ails 
Gauss! 11-ordinate formula is used. The values of @ and t are again taken 
from Table 5; the values of m,(t) are given by [84], in which the y's are read 
from Figure 3 and m,(t) is given in Table 5. 

The pressure distribution can now be obtained from [69]. Graphs of 
p/q are shown in Figure 5. 


* 


la | 1.00 
© Direct Velocity Method 
oO Axial Doublet Distribution Method, Gauss || Ordinate Quadrature 
© Axial Doublet Distribution Method, Gauss 7 Ordinate Quadrature! ——lo.90 
x Kaplan Method 
+ Karmdén Method 


fo) 
h 
(o) 


-0.20 


MO 2 4h a7 LG 65 of aS 02 Ol Be 


Figure 5 - Comparison of Values of p/q 
Obtained by Various Methods 
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ERROR IN DETERMINATION OF p/q 


Let A(p/q), Au, Av, and Am denote errors in p/a, u, v, and m. Then, 
from [69], we have 


at = -2(u Au + v Av) 
from [68], 
A = yl Ale 


and from [67] and [72], except near the stagnation points, 


Au = 240 i sin? 6d6 = 44m 
2 
0 ay 
Hence 
Az = SER) 0 


Bieler 


If now we assume u=1, y' = 0, y® = 4m (Munk's approximation), we obtain 


Thus an error of one percent in the determination of m would introduce an 
error of 0.02 in p/q. 

In the foregoing example the minimum value of p/q was about -0.20. 
Hence an error of one percent in m would have produced an error of ten percent 
in the minimum value of p/q. It was found, in fact, that the results with 
Gauss! 7-ordinate rule deviated from the values of p/q given by the 11-point 
rule by less than 0.003 for the entire body. The /-point rule would have 
sufficed if an accuracy of only 0.003 in p/q were required; see Figure 5. 

If greater accuracy is desired the integrals can be evaluated in the 
forms [71a] and [72a]. If the latter forms are used in conjunction with the 
Gauss quadrature formula the values of x should be chosen identical with the 
t's required by the Gauss formula. This enables the entire calculations, in- 
cluding the iterations and the velocity determinations, to be made arithmeti- 
cally, without resort to graphical operations, so that the method becomes suit- 
able for processing on an automatic-sequence computing machine. In order to 
obtain sufficient accuracy in the integrations and to obtain the velocities 
and pressures at a sufficient number of points along the body a Gauss formula 
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of high order should be used, say n = 16. For this reason the procedure that 
has been illustrated in detail may be less tedious for manual application. 


COMPARISON WITH KARMAN AND KAPLAN METHODS 


In order to compare the accuracy of the Kdrmdn method with the pres- 
ent one, the error function w(x) was computed for a distribution derived by 
the Karman method, employing 14 intervals extending from -0.98£ x £0.98. 
v(x) is graphed in Figure 3. It is seen that the errors are much greater 
than for the present method, especially near the ends of the body. The oscil- 
latory character of (x) is imposed by the condition that the stream function 
Should vanish at the center of each interval. It is conceivable that the 
amplitude of the oscillations in W(x) may remain large even when the number 
of intervals is greatly increased; i.e., the Kaérmdén method may give a poorer 
approximation when the number of source-sink segments is greatly increased. 
The pressure distribution obtained by the Karman method is graphed in Figure 5. 

Kaplan's first method?* was also applied to obtain the pressure dis- 
tribution. Kaplan expresses the potential function @ in the form 


@= FA QL(A) Plu) 


where A and ww are confocal elliptic coordinates, 


P_(u) and QA) are the nth Legendre and associated Legendre 
a polynomials, and the 


An's are coefficients to be determined from a set of 
linear equations which express the condition 
that the given profile is a stream function. 


In the present case it was assumed that g Was expressed in terms of the first 
9 Legendre functions and the An's determined from the conditions that the 
stream function should vanish at 9 prescribed points (including the stagnation 
points) on the body. The resulting pressure distribution is also shown in 
Figure 5. 


SOLUTION BY APPLICATION OF GREEN'S THEOREM 
GENERAL APPLICATION TO PROBLEMS IN POTENTIAL THEORY 


Let @ and w be any two functions harmonic in the region exterior to 
a given body and vanishing at infinity. Then, a consequence of Green's second 
identity?’ is 


ox as - [{u 9? as [85] 
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where the double integrals are taken over the boundary of the body and dn 
denotes an element of the outwardly-directed normal to the surface S. Also 
derivable from Green's formulas is the well-known expression for a potential 
function in terms of its values and the values of its normal derivatives on 
the boundary* 


(a) = qe |[[-2 $2 + 6% Las [86] 


where r is the distance from the element dS on the body to a point Q exterior 
to the body. 

When a distribution of ¢or d¢/dn over the surface of the body is 
given then [85] may be considered as an integral equation of the first kind 
for finding d¢d/dn or ¢ respectively, on the surface. If the integral equation 
can be solved, [86] would then give the value of ¢ at any point Q of the 
region exterior to the body. 


AN INTEGRAL EQUATION FOR AXISYMMETRIC FLOW 


Equation [85] will now be applied to obtain an integral equation for 
axisymmetric flow about a body of revolution. Let y be the ordinate of a 
meridian section of the body and ds an element of arc length along the boundary 
in a meridian plane. Then we may put 


dS = 2my ds [87 ] 


It will be supposed that the body is moving with unit velocity in the negative 
x-direction, which is taken to coincide with the axis of symmetry. The con- 
dition that the body should be a solid boundary for the flow is that the com- 
ponent of the fluid velocity at the body normal to body is the same as the 
component of the velocity of the body normal to itself. This gives the bound- 
ary condition 


‘e = -sin y» [88 ] 


where y is the angle of the tangent to the body with the x-axis. Substitution 
of Equations [87] and [88] into [85] now gives 


12 
[ yo 52 ds = -{ yw sin yds [89] 
0 


where 2P is the perimeter of a meridian section and the arc length s is meas- 


ured from the foremost point of the body. 
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Now let us choose for w an axisymmetric potential function and let 
w(x, y) be the corresponding stream function. Then 


dw _ ay 


y dn ds 


and 


Bs avy a P 
J v9 Gi os = ov] - J ve as 


Also let U be the total velocity along the body when the flow is made steady 
by superposing a stream of unit velocity in the positive x-direction. Then 


ude 
Ws as + cos y 


Furthermore, we have dx = ds cos y, dy = ds sin y. Then [89] may be written 


P P P 
ow =| w(cosy - U)ds = =|] yw dy 
0 0 0 


or 
P 


[90] 


IP P 
j Uwds - | (pdx - ywdy) - ow 
0 0 


0 


But, since w and w are corresponding axisymmetric potential and stream func- 
tions, we have 


Hence wdx - ywdy is an exact differential defining a function Q(x, y) such 
that 


g2 _ y, $e = -yw [91] 
But since also 
y OH . OY 
Via Oox 
we obtain from [91] 
ORO) O29) a2 
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which, by comparison with [4], is seen to be the equation satisfied by the 
Stokes stream function. Conversely, if Q is a function satisfying [4], it can 
readily be verified that the functions w and w defined by [92] are correspond- 
ing axisymmetric potential and stream functions, i.e., that they satisfy Equa- 
tions [3]. Written in terms of 2, [90] now becomes 


oes ds = (@-4 92) [92] 


If we choose for 2 the stream function of a source of unit strength 
situated at an arbitrary point of the axis of symmetry within the body, we 
have, from [10], 


1/2 
TY eer (ecu rere | [93] 
Then 
O62 _ y* 
Ox r3 


and, since y vanishes at both limits, 


(2 - 9§2)[, = 2 


Hence [93] becomes 


jo y* (x) distal [94] 


It is seen that [94] is an integral equation of the first kind in which the 
unknown function is U(x) and the kernel is y*/2r°. 

In contrast with the integral equations for source-sink or doublet 
distributions which can be used to obtain the potential flow about bodies of 
revolution, the integral equation [94] has two important advantages. The 
first is that a solution exists, a desirable condition which is not in general 
the case when a solution is attempted in terms of axial source-sink or doublet 
distributions. The second advantage is that [94] is expressed directly in 
terms of the velocity along the body so that, when U is determined, the pres- 
sure distribution along the body is immediately given by Bernoulli's equation 
[69]. In the case of the aforementioned distributions, on the other hand, it 
would first be necessary to evaluate additional integrals, to obtain the ve- 
locity along the body, before the pressures could be computed. 
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KENNARD'S DERIVATION OF THE INTEGRAL EQUATION 


A simple, physical derivation of the integral equation [94] has been 
given by Dr. E.H. Kennard. This will now be presented. 

Imagine the body replaced by fluid at rest. Let U be the velocity 
on the body. Then the field of flow consists of the superposition of the uni- 
form (unit) flow and the flow due to a vortex sheet of density U. 

Now subtract the uniform flow. There remains the flow due to the 
vortex sheet alone, uniform inside the space originally occupied by the body, 
of unit magnitude. 

A vortex ring of strength Uds produces at an axial point distant z 
from its plane the velocity 


Soe y*Uds 
2(y?+z? ye 2 


where y is the radius of the ring. Let s be the distance of a point on the 
body measured along the generator from the forward end, in a meridian plane. 
The axial and radial coordinates will then be functions x(s), y(s). The ve- 
locity due to the sheet at a point t on the axis will then be 


[SLO as = 
0 


ers 


where r® = [x(s) - t]" + y*(s) and P is the total length of a generator. The 
equivalence of this equation with [94] is evident. 


A FIRST APPROXIMATION 
If we again make use of the polar transformation x - t = y(x) cot 6, 


[94] becomes 


7 U(x) sin? @d@ _ 
o 2 sin[@->(x) ] 


1 [95] 


When x = t, @ = 4. For an elongated body the integrand in [94] peaks sharply 


in the neighborhood of x = t, so that a good approximation is obtained when 
U(x) is replaced by U(t) for the entire range of integration. Also, »(x) will 
be small except near the ends of the body so that the approximation 


sin [6= y(x)]= sin 6 cos y(x) = sin @ cos y(t) 
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will also be introduced. We then obtain from [95] the approximation 
U(t) = cos y(t) [96] 


Just as was done in the case of Munk's approximate doublet distribu- 
tion we can improve upon this approximation in terms of an estimated longi- 
tudinal virtual mass coefficient for the body. For this purpose we will first 
derive a relation between this coefficient and the velocity distribution. 

Let T be the kinetic energy of the fluid when the body is moving 
with unit velocity in the negative x-direction. Then 


dd a 
eT = -p rae dS = amp | y¢@sin yds 
0 


by [88]. Integrating by parts and substituting for d¢é/ds from [93] now gives 
2 dd 12 
QW = -mp | We as ds = 7p | U(x) y*(x)ds - A 
0 0 


where A is the displacement of the body. But also, by definition, 2T = k 4. 
Hence 


P 
A(1 +k) = mp{ U(x) y?(x)ds [97] 
0 


This is the desired relation between k, and U(x). 

Now suppose, as a generalization of [96], that an approximate solu- 
tion of the integral equations [94] is U(x) = C cos y. If this value is sub- 
stituted into [97], we obtain C = 1 + kK, - Hence an improved first approxima- 
tion to U(x) is 


U, (x) 


(1 + k,) cos »(x) [98] 
Equation [98] gives an exact solution for the prolate spheroid. 


SOLUTION OF INTEGRAL EQUATION BY ITERATION 


In order to solve [94] by means of the iteration formula treated 
in Reference 17, it would be necessary to work with the iterated kernel of 
this integral equation. Since this would entail considerable computational 
labor it is proposed to try a similar iteration formula, but employirg the 
original kernel: 


(t) = U(t) + cos y(t) |1 - pe U (x)as] [99] 
0 r 


Une n 
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where r= = (x - t)? + y*(x) and x = x(s). 
Here also it is convenient to express the iterations in terms of 


error functions E(t) defined by 


PIAUIN (xe) mays (9) 
us = n 
E(t) = 1 j aaa ds [00] 


or, from [99], 


E(t) cos y(t) = U,,(t) - u(t) [101 ] 
Hence 
n 
Unsq(t) = U,(t) + cos »(t) & Flt) [102] 
Also, from [99], 
2, E (x) y* (x) 
Bq (t) = E(t) tae aaa [103] 


where X,» X, are the nose and tail abscissae. Thus, to obtain Una ft), we 
first obtain E(t) from U, (t) in [100], then 5 i , Be from [103], and 
finally Une 't) from [102]. 
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NUMERICAL EVALUATION OF INTEGRALS 


In applying Equations [100] and [103] it will frequently be neces- 
sary to evaluate integrals of the form 


(See) dx, where r® = (t-x)® + y?(x) 
Zo r 


This form, however, is unsuited for numerical quadrature for elongated bodies, 
since y*(x) peaks sharply in the neighborhood of x = t. Here, as in the case 
of the integrals for the doublet distribution, two procedures are available 
for avoiding this difficulty. The first employs the polar transformation [70], 
involves several graphical operations, but in general transforms the integrand 
into a slowly varying function so that the integral can be evaluated by a 
quadrature formula using relatively few ordinates. The second retains the 
original variables and eliminates the peak by subtracting from the integrand 
an integrable function which behaves very much like the original integrand in 
the neighborhood of the peak. The numerical evaluation of the resulting inte- 
gral on the second method requires a quadrature formula with more ordinates 


RAR 


4y 


than the first in order to obtain the same accuracy, but, since all graphical 
operations are eliminated, the second method is suitable for processing on an 
automatic-sequence calculating machine. 

The result of the polar transformation has effectively been given 


n [95]. We have 


“1 E(x) = x) E(x) sin? @ ai pe ) 4 
J x= | sin[6 ->(x) [104] 


0 r° 


where 


x - t = y(x) cot 6 [70] 


It is desired to evaluate this integral for a series of values of t. In gen- 
eral this can be done with sufficient accuracy by means of the Gauss 7- (or 
11-) ordinate quadrature formulas. This gives 7 (or 11) values of 6 at which 
the integrand needs to be determined for a given t. The value of x occurring 
in the integrand is determined implicitly, for given values of t and 6, by the 
polar transformation [70]. In practice the 7 (or 11) x's can be obtained 
graphically from the intersections with a graph of the given profile of the / 
(or 11) rays from the point x = t on the axis at the angles required by the 
Gauss quadrature formula. If greater accuracy is desired, these graphically 
determined values of x can be corrected by means of the formula 


t-x_ + y(x_) cot 6 
xe bo, See ae [105 ] 
& I-y! (x,) cot @ 


in which x_ is the graphically determined value and y!' denotes the derivative 


of y with respect to x. 
Now let us derive an alternate, completely arithmetical procedure 


for evaluating the integrals. Put 


k' (x, t) = Fale t) 
[Gxt Pe eles) RY 
where y> = f(x) is the equation of the given profile and y* = g(x, t) is 


the equation of the prolate spheroid whose ends coincide with the ends of the 
given body, and which intersects the given body at x=t. i-.e., 
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(x-X 9) (x, -x) 
g(x, t) = f(t) (t-x,)(x,-t) [106] 


The length-diameter ratio A of the spheroid is given by 


(t-x 4) (x, -t) 


ee ea ee [107] 


whence the longitudinal virtual mass coefficient k(t) ean be obtained from 
[26]. 

Since U(x) = (1 + k,) cos y(x) is an exact solution of [94] for the 
prolate spheroid, we have 


Jets, t) ax = Te TE] [108] 


We now obtain, from [98], [100], and [108] 


1+k 
1 


Pree [109] 
1+k, (t) 


[ tetx, t) - k'(x, t)]dx - 
i) 


Also [103] may be written in the form 


n+1 


By i(t) = Balt) - 5] k(x, (E(x) - Eye) Jax - B,(t)[ k(x, tax 
0 0 


But from [98] and [100], 


1-E, (t) 


Za 
Nese, ie )lebe = 2 === 
i: THe 


Hence we obtain 


[etx, t)[E (x) - By(t)jax [110] 


pean ameTegm ata ad 
0 


ILLUSTRATIVE EXAMPLE 


The present method will now be applied to the same profile [78] as 
before. By way of contrast with the semi-graphical procedures previously used, 
a completely arithmetical procedure will be employed. 
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The velocity U(t) will be determined at the 16 points along the body 
whose abscissae are ty = by the Gaussian values for the 16-point quadrature 
rule, Table 1. Since the body is symmetrical fore and aft, it is necessary to 
determine the velocity at only half of these points. Values of y(x), cos y(x) 
and k, (t) for these points are given in Table 6. 

In order to apply the Gauss 16-ordinate rule it is necessary to eval- 
uate the integrands in [109] and [110] at the 16 Gaussian abscissae x, = & 
for each of the 8 values of t,- Thus, there are 16 x 8 = 128 values of k(x, t) 
and of k'(x, t) to be determined. The matrices Ki = R k(x, t,) 
Lae = Rik (x5, t,) where the R,'s are the Gauss weighting factors, are given 
in Tables 7 and 8, and applied to evaluate E(t) from [109]. E,, Ej, and oy 
are then obtained from [110]. u(t) is then given by [102] and then p/q by 
[69], in the form p/q = 1 - De The arrangement of the calculations and the 


and 


results are given in Table 9. The graph of p/a is included in Figure 5. 


TABLE 6 


Values of y, cos y, and k, (x) for Application of 
Gauss 16-Point Quadrature Formula 


Psa [rr |v Prove 


-0.9894.009]0.0408548] 1 .8965483) 1.0856 QO. Se: 


9445750] .0903198]0. 7464764) 0.6412 093389 
.8656312| .1324422) .3917981 .088359 
27554044) .1642411] .2099651 .081862 


.6178762} .1848527 
4580168] .1955501| .0393076 
.2816036| .1993706| .0089607 


-0.0950125 |0.1999919]0.0003431}0.0003431 


. 1020867 .074689 
.067885 
.062506 


0.059509 


44 


TABLE 7 
f(x,) 


Matrix of Values* kK,. = R, = 
ji J p 2 3/2 
[(«; t,) + t(x,)] 


0715 
- 


.0098 

.00505 

.00341 

.00208 

.00121 

.00062 

00027 ‘4 
1 Sie OOO OMe .00009 3 2 
1610.00001}0.00001}0.00001| 0.00001 |0.00001| 0.00001 | 0.00002 | 0.00004 


TABLE 8 
Ban ea )) 
Matrix of Values** in = R ——__— _=_____ 
[Cx v7 t,) +e(x5, t,)] 


0.00003 
.00044 


.05223| .06265| .08976] .16366} .40559| .86501| .34633 


03061} .03502| .04556] .07041} .13995! .39080] .91586 
01805] .01998] .02432] .03352| .05520} .11924| .37496 
01058 .011 01329] .01094| .02457] .0 Nf 10217 

| .03396 
11| .00334] .00351| .00385 oon] .00561| .00787] .01303 


3 

ie .22699| .4 : 33504) 07174 I. 
Be 11637] .19118] .418621 .80927] .34255] .07641 

9 


10] .00607| .00645} .00725] -..00876] .01165] .0178 


12 .00178 .0021 .00262| .00346] .0052 
00081] . 


15| .00007 
1610.00001| 0.00001] 0.00001] 0.00001] 0.00001] 0.00001] 0.00001] 0.00002 


*For i>8useK.. =K : ; 
Ji dP a) 5. dbt(at 


**For i > 8 use K'. = ail ‘ 
ae lag IRL 
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TABLE 9 
Calculations for E(t) and U(t) 


Assume k; = 0.06: Put Boy = Rk 7 Si Simi 15 (ze sa (e,)) =D 


ji fil dj ag a” ial} nj 
(a)  S OnSRSONs aes p & Og —— A 6. 
a ab . 3 v= Vo Tera (Gay) T . 2) 


ae eee re) 0.66458 
28oee 


.00118 .00566 
.00078 .00706 
.00048 00802 
.00028 .00851 
-00851 
.00802 
.00706 


00 3 
0 


0 Bilin : : 7 

: -00171 -0520 

-00062} .00078 ~00016 -03333 
"4 -00027} .00029 =0. ae .01799 

-00007} .00007 -0.00730 
12 0.00001 |0.00001 0} 


Us (x, ) = 0.5450, 


lt+tk 


(2) 02S OLMIED C287? <2) as Oly T—a 


0.96947 


1 }0.20313]0.20017| +0.00296 
. 68926} .6892 -0.00001 
+0. aos 


.21382|_. i 
‘04486 .06265 :01779 
.02187| .03502 .0131 
-01140|_ .01998 .008 0010 02 - 00029 
.00639| .01143 0050 09357 00060 02563 00016 
.00379| .00645 . 00266 .0 .00031 02312 . 00009 
.00228] .00350 .00122 064.83 .00015 .01875 . 00004 


.00131| .00178 0004 OU -00006 «01 344 . 00002 
02603 5 


.000 067 . 00081 - 0001 
$0028 soueee -0. eggiee 
.00008] .00008 
0.00001 0.00001 
Be 


46 


TABLE 9 (Continued) 


l+k 


(en) Oe pees aeceihs Oo = Oe at Taco 


= 0.97391 


+0.00206 } +0.0003 : +0.00013 : +0. 0000 
+0.00381 0.0106 +0.0031 A +0.00096 : +0.0002 
-0.00001 
+0.00240 


5} .00009| . +0.00001 +0. 00327 
16]0.00001 {0.00001 0) +0.00634 


= -0.06312 -0. ras = -0. OG 
3 = 7005951 = +0.01546 = +0.00427 


E 
sb) 
ra +0.11275 


l+k 


0.00338] 0.00271 +0.00067 +0.03333 +0.00011 +0.00004 +0. 00403 +0.00001 
04802 .00712 02603 007 4 : 00044 00245 .00014 
00465 01534 +0.00483 +0. +0.00152 +0.00150 +0. 00049 
0) 0) ) 0) 0) 
00068 One =0. 00783 =05 =0.00230 =0. 00165 -0. 00069 
.02019 038 00557 : poole 00303 
.01697 .05666 .00303 ; .00081 “oot 
01045 06754 ; 6 0176 00041 -OOKY, 
-00557 20075 .00077 : 
.00260 ; -00035 
00099 , .00014 
-0.00026 ; -0.00004 -0. 
0 0 0 
+0.00004 onsen +0. pod +0. +0.00150 
+0. poco! . 02603 00245 
0.00001 {0.00001 +0.03333 +0. ae +0. 00403 


-0. Ree -0.01290 = -0. aa -0. ere 
= +0.04u17 = +0.01079 = +0.00277 = +0.00071 


= 0.09827 U(x, ) = 1.0946; + = -0.1981 


See 
1 
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TABLE 9 (Continued) 


l+k 


Gatsepi; 


(e) x, = -0-617876; cos 7 = 0.9948; 


+0.00026 +0.0000 j +0.00001 +0. 00568 0 
ee Uy 0005 ; .00018 0047 0 +0.0000 
011 : 00257 i 000 00315 0002: 

00616 4 4 +0.00641 ; aE +0.0016 +0.00056 


0 

-0.00138 -0.0005 

00234 0002 

; 0028 00012 
: : 700023 283 

"03792 ; 00009 

-0.02006 0.00010 ; -0.00003 

0 0) 0) 0 

+0.00015 0187 : : +0. 00001 +0. 00165 

00011 ; ; +0. goat 00315 

+0. Gyo 07 344. 00410 

+0.01651 +0.00568 


1 
E 
3 
i 
: 
: 
10 
1 
12 
13 
14 
15 
16 


U_(x_) = 1.0864 ; 2= -0.1803 
5.5 q 


1 + ky 
f X_, = -0.458017; cos y = 0.9992; —————— = 0.99262 
(2) 36 [1 +m (x@)] 
7 
EB. =p me (i. 13 }) E_ -E K -E__) 

2j 26 je 23 26 33 36 age eal <3 

+0.00012 4 +0.00002 +0.02182 +0.00001 +0.00706 0 
.00146 ; Cee aoe) . 00008 00548 +0.00002 
01548 -00033 00453 .00010 
01081 ata 0030 00026 
+0. 00537 +0.001 +0. 0013 +0. 00048 

0) 0) 0 0 
=. -0.00168 -0.00096 -0.00037 
0068! oe | .007 4 00015 
00688 0002 00145 700005 
-0.00437 -0.00007 -0.00096 -0.00002 

0) 0 0 0 
+0.00008 +0.00531 +0.00002 +0.00138 +0.00001 
a; ; .00007 sae a 00002 700305 +0. 00001 

i ; 00004 : +0.00001 00453 0) 

j 3 +0.00001 : 1875 0) 00548 0) 

i PO 0 +0.00706 0) 

= +0.00379 f = +0. = +0.00052 i = +0.00029 
E_ = +0.00537 E_ = -0.00002 E_ = -0.00026 E = -0.00016 
16 26 36 46 


+ E 
eas = = ; B=. 3 
SSE = ORGY Us\(ocs) 1 06H e025 
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TABLE 9 (Continued) 


l+k 


(g) 7 = -0.281604; cos y = 1.0000; ea 


= 0.99764 


+0.00000 
-00001 
. 00004 
-00011 


+0 .00000 
00006 
00006 


“00002 
+00001 
+0.00000 


7) = 


(h)  x,= -0.095013; cos y = 1.0000; 


k! 


0.00003 0.00003 0.00000 0.00851 0.00000 

00044 00037 .00002 00693 00001 

ony .00008 00 fe .00002 

00 .00019 00448 0000 

Scone 00283 0000! 

01079 00145 .00013 

00507 00049 .00917 

00000 .00000 .00000 

00000 00000 .00000 

00420 00049 00005 

00456 00145 .0000' 

00321 0028 .0000 

0051 00339 00179 0044 00002 
.00183| .00107 00076 00548 00001 
15 OE .00024 00021 : .00693 .00000 
16|0.00004] 0.00002 0.00002 0.00000 0.00851 0.00000 

J = +0.0432 J = +0.01219 f = +0.00063 
[Bs = -0.0233 -0.00690 -0.00038 
k 
i 0.03456 
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SUMMARY 


Two new methods for computing the steady, irrotational, axisymmetric 
flow of a perfect, incompressible fluid about a body of revolution are 
presented. 

In the first method a continuous, axial distribution of doublets 
which generates the prescribed body in a uniform stream is sought as a solu- 
tion of the integral equation 


where r is the distance from a point (t, 0) on the axis to a point (x, y) on 
the body, r= = (x - t)® + y=(x). 

A method of determining the end points of the distribution and the 
values of the distribution at the end points is given. If the equation of the 
body profile, with the origin of coordinates at one end, is 


ai 2 3 
uP (8) SEL Ok te BS ar ELIS HP Goc 


a very good approximation for the distribution limit a at that end, when the 
coefficients a,, a5, .-- are small, is given by 


a 
1 1 
ae y+ a, +> Vala, 
if a, 2 0), 162 a, is negative, the term containing it is neglected. The cor- 
responding value of the doublet strength at this point is 


Inia a8 uD } o§ a Vaia, 
Formulas and tables for determining a and m(a), which may be used when the 
above procedure is insufficiently accurate, are also given. The values a, b, 
m. = m(a), m, = m(b), f= y-(a) and fe y-(b) are then used to obtain the 


a 
approximate solution of the integral equation 


i fy (Date a ey SEE ames — 
mei X) =" Cl Ven oa pend Sehr! pacacn Wid pen eat a DES 
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where 


1+k, (71 

; m 2 | eax - 3(b-a) (m, +m, ) 

Tea Oe Re SE ee 
i yrdx - 5(b-a) (f,+f,) 


and k, is the longitudinal virtual mass coefficient for the body. 
This approximation is used to obtain a sequence of successive approx- 
imations by means of the iteration formula 


B i, ((te)) 
my (x) = (x) +3 v(x)/p - [ + at] 


When a doublet distribution has been assumed, the velocity components at a 
point (x, y) in a meridian plane are 


and the pressure is given by 


where q is the stagnation pressure. 

The iterations are most conveniently performed in terms of the dif- 
ferences between successive approximations to m(x), which also furnish, at 
each iteration, a geometric measure of the accuracy of an approximation. 
Simpler forms for the velocity components at the surface of the body are given 
in terms of this difference or error function. 

Gauss! quadrature formulas are recommended for the numerical eval- 
uation of the integrals. Two methods of carrying out the iterations are 
given. The first employs a polar transformation and a graphical operation be- 
tween successive iterations; the second is completely arithmetical and is 
suitable for processing on an automatic-sequence computing machine. All of 
these procedures are illustrated in detail by an example, in which the semi- 
graphical method is employed. The accuracy of the method is analyzed; the re- 
sults are compared with those obtained by the methods of Karman and Kaplan. 
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In the second method the velocity U(x) on the surface of the given 
body is given directly as the solution of the integral equation 


[Sec ieee 
0 ae 


where s is the arc length along the profile, 
x is equal to x(s), and 
2P is the perimeter of a meridian section. 


An approximate solution to this integral equation is 


U, (x) = (7 + k,) cos y(x) 


where kK is the longitudinal virtual mass coefficient and y = arctan — 


U, (x) is used to obtain a sequence of successive approximations by means of 


the iteration formula 


U 


nay (t) = U,(t) + cos y(t) | “(2 U,(x)¢s 


Here, also, the iterations are most conveniently carried out in terms of the 
differences between successive approximations to U(x) which also furnish a 
measure of the error in the integral equation. Two methods of carrying out 
the iterations are again available, of which one is semi-graphical, the other 
completely arithmetical. The latter technique is employed on the same example 
as was used to illustrate the first method. 
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APPENDIX 
END POINTS OF A DISTRIBUTION 


An approximate method for determining the end points of a distribu- 
tion and its trends at the ends will now be described. Let y* = f(x) be the 
equation of the given profile extending from x = 0 to x = 1; let m(x) be the 
corresponding doublet distribution, extending from x = a to x = b. It will be 
assumed that 0 < a << b < 1 and that a is near 0, b is near 1. Then m(x) is 
given by the integral equation 


(eit) 


Various conditions on m(x) may now be obtained by differentiating 
[111] repeatedly with respect to x. We get 


f B@ltex - 26 + e1(x)] at = 0 [112] 
a r? 
: 5 1 1 " , 
Ji m(e)[-S (2x - 2b + 8 2 Lieve )Jat = 0 [113] 
[ince 22 (2xa2t- tun) oe 1) (24£") (2x-2t+f!) n Aas = 0 [114] 
a Uy? er! r> 


When x = 0, r = t and, writing f(x) as a Taylor expansion 


(x) 


2 3 
Anos ap BLO ar BL Se SP oo c [115] 


then also f'(0) = a,, f¥(O)) s 2a,, oY (O) Ss 6a. Now, setting x = 0 in Equa- 


tions [111] and [113], we obtain 


b 
m(t ) ull 
it = at = 4 [116] 
6 
f mela, SUE atH= nO 117] 
b 
I a Se, - 20a,t + 4(4 - a, )t?|at = 0 [118] 


D2 


9 


. me [58s - 210a,t + 60a, (6-a,)t? + 4O(3a,-H)t? + ala,t*|at =o [119] 


Also assume that m(x) may be expressed as a power series 
= 2 
m(x) = ¢, + c,x + Ce ap boc [120] 


Then Equation [116] gives 


HL -L) 4 eft-4)4 0, ms2e 


or, neglecting 1/b° in comparison with 1/a* and setting b = 1 in comparison 
walide Wily 


2 al ee 
tp 2c a(1-a) + 2¢,a" log A oad = [121] 
Similarly, Equations [117], [118], and [119] give, approximately 


c, (3a, -8a) + 4e ala, -3a) + 6c,a"(a, -4atta*) = 0 [122 ] 


2c, | 5a, -2a,a+6(H-a, )a®| + Ye a/3a -15a,ar4(4-a, )a*| 


[123] 
+0,8°| 15a -80a,a+2u(4-a, Ja? | = 0 
1 -240a,a+80a, (6-a, )a?+6l (3a, -4)a9+48a,a*| 
30, |35a, 240a,a4 Oa, ( -a,)a + (Ja, Ja“+48a,a 
+2lte, |5a,a-35a,a2+12a, (6-a, )a°+10(3a,-4)a*+8a,a>| [124] 


the, [35a,a?-252a,a°+90a, (6-a, )a*+80(3a,-4)a°+72a,a°] = 0 


Equations [121] through [124] are sufficient in number to determine the un- 


knowns a, c Since the latter three equations are linear and homo- 


Cc Cc 
Oe stole 
geneous in Cro Gap and c,, a can be determined from the condition that the de- 
terminant of their coefficients must vanish. In this way the following equa- 


a 
tion of the 7th degree ina => was obtained: 
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a(a - 4)?(5a* - 830° + 28807 - 3680+ 128) - 96a or( oe - 4) 
+ 4a ala - 4)(5307 - 148a + 128) + 1152a, a, (2c - 3) 
[125] 
+ 72a, (a - 4)? (50° - 250° + 40a - 16) + 48a,a,a(3e - 8) 


2 
= Asc e C!S M\(Gee— NS slo) = Vis2e eo oY) = © 


Corresponding to a solution a of [(125], Cy» Cy» and c, can be obtained from 
Equations [121], [122], and [123]. The solution of the latter equations gives 


¢,D = -4a?| 30° Sr tO SO ey cp Beles = Gear 10 = Ya. )] [126] 


e,D a|150° = ser ts Glee S oehi ae Glee ap UtRW(Ger = ais At & 6a,)) [127] 


e,D = -4[(a - 4)2(a - 1) + Hag] [128] 
where 


D = 2(9a° - 94a? + 272a - 192) + 8| (a - 4)?¥(e - 1) + Ya, | 1og a + 96a, 
-2a(15a° - 264a° + Guta - 768) - 384aa, - 96a*(5a° - 2ha + 2k) 


+ 576a7a, [129] 
The initial doublet strength at x = a is 
= 2 
m(a) = ¢, Uaaeene CLES Se lows 


or, from Equations [126] through [129], 


m(a) = -2-[(a - 4)(a? - 120 + 16) + 4Bala - 4)(a - 2) + 16a, - 96aa,] [130] 

' Equations [125] through [130] determine the end points of the distri- 
bution and its initial trends. In general, Equation [125] will have more than 
one real root. In this case the initial trends corresponding to each of the 
roots should be examined, and that root chosen which appears to give the 
"simplest" trend. 


DD) 


The equations can be solved explicitly in the case of a very elon- 
gated body for which a,, a,, a,, ... in [115] are all very small. First let 
us suppose that they are so small that all the terms in [125] containing them 
are negligible, so that the first product term alone may be equated to zero, 
abo San 


ala - 4)?(5a* - 830° + 28807 - 3684+ 128) = 0 [131] 


whose real roots are a = 0, 0.547, 4.0, 4.0, and 12.429. 

Let us consider the solution w = 4; i.e., a= ie Since the radius 
of curvature at x = 0 is avV/en this solution is seen to be in accord with 
Kaplan's assumption for the end points of the distribution. Furthermore, sub- 
stituting w = 4 into Equations [129] and [130], we obtain, to the same order 


of approximation, 


a a, 
D = 64, Gy = Tas Os img Cie 0) 
whence 
2 
Ea Sa 
m(x) = Rare ae Te 8 m(a) = 0 [132] 


In order to obtain a second approximation it will be assumed that 
not only a,, 43, 4, ... but also (@ - 4) are small to the first order. Then, 
neglecting terms of third and higher order, Equation [125] becomes 


-3072(a - 4)? + 6144a (oe - 4) - 3072a, + 768a,a, = 0 [133] 
whence 
a=4+a t 1 yaa, [134] 
provided 
a 2 0 


Corresponding to this value of a we obtain from Equations [126] through [129], 


to the same order of approximation, 


*The smaller of these two roots has given the preferred solution in all cases tried thus far. 
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2 


ay i 
m(x) = c(-— +ax+a x + Lg 


where [135] 


ib) 
fo 


cad + +5 10) 


and 


m(a) == Ca" yaa, [136] 


The expression for m(x) in [135] may also be written as 


2 


n(x) = (-,2 + ¥) [1354] 


When a, < 0 the solution for a in [134] indicates that there would 
be no real roots near a = 4. In this case a graph of the complete polynomial 
in [125] should be examined either for the possibility that more complete cal- 
culations would show that there are real roots near a = 4 nevertheless, or 
that the maximum value of the complete polynomial in the neighborhood of o = 4 
is so nearly zero, that the value of @ corresponding to this maximum may be 
taken as an approximate solution. On this assumption, the second order analy- 


sis would give 
Gye ob AP Bop a, < 0) [137] 


Since a, does not occur explicitly in Equations [135], it is seen that they 
would also be obtained, to the same order of approximation, if the value of 
o in [137] were substituted into Equations [126] through [129]. 

If it is determined that not even an approximate solution can be 
assumed near ow = 4 it would be necessary to consider solutions in the neighbor- 
hood of the other roots of Equation [131]. 

In order to facilitate the computations for graphing the polynomial 
in [125], the functions A(a), Bla), ... H(a), where 


a1 


A(a) = ala - 4)?(5a* - 830% + 288a7 - 368a + 128) 
Bla) = 72(a@ - 4)?(50° - 250° + 40a - 16) 
C(a) = Yala - 4) (5307 - 148 + 128) 


D(a) = -288(a@ - 4)(5a® - 160 + 16) 
[138] 
E(a) = -96a(3a - 4) 
F(a) = 1152(2e - 3) 
G(a) = 48a(3a - 8) 
H(a) = -1152(@ - 3) 


have been tabulated in Table 10. In terms of these functions, Equation [125] 
becomes 


2 2 2 
A +a,B+a,C + a,a,D + aj,E + a.a,F + a,a,G + a,a,H = 0 [139] 


It is of interest to compare the approximate value for oa from Equa- 
tion [134] with the exact value for the prolate spheroid y* = sab: =e) ein 
this case we have 


1 
Seer Te2en nee 3 


and the exact value of a is 


1 1 1 
mas 2y/1 a ape Rca Pec a 
NE a un* 


But when the length-diameter ratio A is large, Equation [134] gives the ap- 
proximate value a = 4 = sa which is seen to consist of the first two terms of 
the series expansion of the exact value of a. Table 11 shows that the approx- 
imate formula gives excellent agreement with the exact values even for very 
thick sections. Both the exact and the approximate formulas give m(a) = 0. 
Thus the present approximate methods work very well for the prolate spheroid. 
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TABLE 10 


Functions for Determining Limits of Doublet Distributions 


[a [xe [ote [at [we rf ro fot Pn] 


(a) 
3456.0 


-13409 S bs 16230 B. -3225 
.6| 14227 65 .3|-2995 
Al 12414 89 .3|-2764 
.1] 10782 107.5] -2534 


0 

.6 

3 

8 

4 
9324 -2304.0 2880 
8029 -2073.6 
6890 -1843 32 264 
5898 -1612.8 253 
5044 -1382 .4 2419 
4320 -1152.0 2304 
3716 -921.6 2188 
3225 .2 2073 
2838 a} 1958 
2545 4 1843 


MEDMO NERHDO NEDHMO NEDHHDO NERHHMO 


0 0 
0 7 a2 
5 5 3B 
si an .2 
5 9) off 
4 33} 8 0 
4A 6.6 .6 4 
ou oll 4 4 
5 4 6 32 
8 1 al 03 
0 0 0 0 
.2 03 3 6 9 
5 of J 6 ol 
9 33} 4 2 5 
9 23) 03} 9 9 
9 33 .8} 2340.0] -72.0 1728 
ol 9 .9] 2211.8] -122.9 4 
5 5 2 2152.8] -179.5 8 
32 at : 2154.2] -241.9 .2 
8 4 -6| 2207.5] -310.1 6 
.0 0 .O| 2304.0} -384.0 .0 
8 .0 .8] 2435.0) -463.7 4 
=) mil -3] 2592.0} -549.1 8 
ae} 0 -1| 2766.2) -640.3 .2 
.2 ot .3} 2949.1] -737.3 6 
23 8 .8| 3132.0] -840.0 .0 
0 29 2 3306.2] -948.5 4 
‘f 5 : 3463 .2]-1062.7 8 
3 9 4) 3594.2)-1182.7 o2 
29 9 .2| 3690.7)-1308.5 .6 
0 0 .0] 3744.0}-1440.0] 3456.0 
8 .2 4] 3745.4)-1577.3] 3686.4 a2 
ol .2 .5| 3686.4)-1720.3] 3916.8 4 
6 8 .0| 3558.2]-1869.1] 4147.2 6 
9 3 .8] 3352.3|-2023.7| 4377.6 8 
03 03 .8| 3060.0]-2184.0] 4608.0 0 
o® 5 .8] 2672.6]-2350.1] 4838.4 ip 
8 of .5| 2181.6]-2521.9] 5068.8 4 
8 2 .0 Taubes -2699.5| 5299.2 6 
0 6 8 53 .9|-2882 .9] 5529.6 8 
-3072.0 0 0 
5 9 .2| -3266.9 4 2 
5 8 2) -3467.5 8 4 
5 33} -8] -3673.9 22 6 
oD 3 -5|-3886.1 6 8 
8 03} -0]-4104.0 0 0 
.6 B2 .0}-4327.7 4 12 
2 otf -0}-4557.1 8 4 
4 é 8 03 2 6 
5 9 a5 6 8 
0 


offo NM ENaON NW EN 


10 | -2790720 -615168 
11 | -3417260 -897120 
12 | -1966080 -1253376 
13 | 4706910 3351348] -1692576 


-2223360 
-2854368 


22052800 4728640 


A+a,B+a,C+a.a,D + a,E + ajaF + a,a,G + ava 


A = ala - 4)?(5a* - 830° + 288° - 3680 + 128) E = -96e(3a - 44) 
B = +72(a - 4)? (5a° - 250° + 40a - 16) F = 1152(2a - 3) 
C = Hala - 4) (530 - 1480 + 128) G = 48(3a - 8)a 
D = -288(a - 4)(5a* - 16a + 16) He se 6 3) 
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TABLE 11 
Comparison of Exact and Computed Values 


a 
of a = aa for a Prolate Spheroid 


A 2 
Exact o Dol D2 |B 


Approximate a@| 3.750 


(eo) 
On 
WW 


\N 
No} 
o 
re 
\N 
No} 
= 
DO 
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